Universal shocks in random matrix theory 
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We link the appearance of universal kernels in random matrix ensembles to the phenomenon of shock forma- 
tion in some fluid dynamical equations. Such equations are derived from Dyson's random walks after a proper 
rescaling of the time. In the case of the Gaussian Unitary Ensemble, on which we focus in this letter, we show 
that the orthogonal polynomials, and their Cauchy transforms, evolve according to a viscid Burgers equation 
with an effective "spectral viscosity" v, = 1/2A', where N is the size of the matrices. We relate the edge of 
the spectrum of eigenvalues to the shock that naturally appears in the Burgers equation for appropriate initial 
conditions, thereby obtaining a new perspective on universality. 

PACS numbers: 



In a seminal paper fT], Dyson showed that the distribution 
of eigenvalues of a random matrix could be interpreted as the 
result of a random walk performed independently by each of 
the matrix elements. The resulting distribution yields the so 
called "Coulomb gas" picture, with the eigenvalues identified 
to charged point particles repelling each other according to 
Coulomb law. For matrices of large sizes, this correctly de- 
scribes the bulk properties of the spectrum |2]. In his original 
work Dyson introduced a restoring force preventing the eigen- 
values to spread for ever as time goes. This is what allowed 
him to find a stationary solution corresponding to the random 
ensemble considered with a chosen variance (related to the 
restoring force). In this note, we point out that a richer struc- 
ture emerges if one ignores the restoring force and performs a 
rescaling of the time of the random walk. The random walk 
can then be described by an equation of fluid dynamics, the 
viscid Burgers equation. In this picture, that we may refer 
to as "Dyson fluid", the edge of the spectrum appears as the 
precursor of a shock wave, and its universal properties follow 
from a simple analysis of the Burgers equation that was devel- 
oped in other contexts 1 3] . This enables us to recover familiar 
results of random matrix theory in a simple way, and get a new 
perspective on the issue of universality. 

In this letter we shall consider only the Gaussian Unitary 
Ensemble, although we believe that many of our results can 
be extended to other ensembles. Thus, we consider N x N 
hermitian matrices H with complex entries. We assume that 
these matrices evolve in time according to the following ran- 
dom walk: in the time step At the matrix elements change ac- 
cording to Hij ^ Hij + 6Hij, with {6Hij) = 0, and {(6Hijf) = 
(l+6ij)At. That is , we assume that at each time step, the incre- 
ment of the matrix elements follows from a Gaussian distribu- 
tion with a variance proportional to Af. The initial condition 
on the random walk is that at time f = 0, all the matrix ele- 
ments vanish. Alternatively, let x, denote the eigenvalues of 
H. The previous random walk translates into a corresponding 
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random walk of the eigenvalues, with the following character- 
istics lilj 



{6xi} = E(xi) Af, 
where the "Coulomb force" 



{(dxif) = Af, 



(1) 



(2) 



originates in the Jacobian of the transformation from the ma- 
trix elements to the eigenvalues, A = n(<>(^/ ~ ■^/)^- Ths 
probability for finding the set of eigenvalues Xi, - ■ ■ at time 
t obeys the Smoluchowski-Fokker-Planck equation 



dP 

'dt 



_ 1 y _ y 



■ (E(xdP) , 



2 ^ dx^ ^ dxi 
whose solution reads 

P(xi, - ■ ■,XN,t) ^ cY\(xi - Xjf e" 2' ^ 



(3) 



(4) 



with C a (time-dependent) normalization constant. 

The average density of eigenvalues, p(x), may be obtained 
from P by integrating over - 1 variables. Specifically: 



/N N 
Y\ dxk P{xi, ■ ■ ■ , Xa,, f) ^ S(x - xi) 
k=\ 1=1 



(5) 



with normalization J dxp(x) - N. Similarly we define the 
"two-particle" density pix, y) ^ {'Zti ~ ^')'^0' ~ ^i))' 

with J dxdy p(x,y) - N(N - 1). These various densities obey 
an infinite hierarchy of equations obtained form Eq. ([3]) for P. 
Thus, the equation relating the one and two particle densities 
reads 



dp(x, t) _ 1 d^pix, t) 



OA J 



dy 



p(x,y,t) 



dt 2 dA^ OAJ - x-y 
where ^ denotes the principal value of the integral. 



(6) 
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In the large limit, this equation becomes a closed equa- 
tion for the one particle density. To show that, we set 



p(x,y) = p{x)p(y) + pco„{x,y), 



(7) 



where pconix,y) is the connected part of the two-point den- 
sity. Then we change the normalization of the single particle 
density, defining 



p(x) = Np{x), 



(8) 



and similarly pix,y) = A^(A^ - 1 
rescale the time so that t - Nt 



'(x,y). At the same time, we 
One then obtains 



dpjx) 
dT 



dy 



P(y) ^ 1 dMx) 

X — y 2N dx^ 



f 



dy 



Pcon(x,y) 

x-y 



(9) 



In the large limit, the right hand side is negligible, leaving 
as announced a closed equation for p{x, r). This equation can 
be further transformed into an equation for the resolvent 



G(z,t) = (^Tr 



1 



z- H(t) 



dy 



z-y ' 



(10) 



whose imaginary part for z - x- ie, and x real, yields the aver- 
age spectral density p(x), and the real part is the Hilbert trans- 
form of p, 'Hp{x) - j- dypiy)/(x - y). By taking the Hilbert 
transform of Eq. (|9]l, keeping only the dominant term in the 
large limit and using well known properties of the Hilbert 
transform, one gets the following equation for G{z, t) 



drGiz, T) + G(z,T)dfi(z, t) = 0. 



(11) 



This is the inviscid Burgers equation |6]. Note that the Lapla- 
cian, that naturally appears in the description of diffusive pro- 
cesses, has disappeared in the large A^ limit, as already indi- 
cated after Eq. (|9j. The non-linear term plays the role of an 
effective mean field representing the mutual repulsion of the 
evolving eigenvalues. We shall return later to the role of dif- 
fusion, and focus now on the solution of Eq. (fTTT l. 

This can be obtained by using the method of (complex) 
characteristics |9], with the characteristics determined by the 
implicit equation z - ^ + tGo(^), where Go(z) - G(z, t - 
0) - l/z- Assuming the solution ^(z, t) to be known, the 
Burgers equation can be solved parameterically as G(z, t) - 
Gq(^{z, t)) = Go(z - tG{z, t)). The solution of this equation 
that is analytic in the lower half plane is 



G(z,t)^^(z- ^|z^-4T), 
2t 



(12) 



whose imaginary part yields the familiar Wigner's semicircle 
for the average density of eigenvalues. In the fluid dynam- 
ical picture suggested by the Burgers equation, the edge of 
the spectrum corresponds to a singularity that is associated 
with the precursor of a shock wave, sometimes referred to as 
a "pre-shock" |3]. This singularity occurs when the mapping 
between z and ^ ceases to be one-to-one, a condition required 
for the validity of the method of characteristics. This takes 
place when dz/d^ = 0=1-1- tG[^(^c), defining ^e(T"). Since 
G;,(^,) = -l/f-, ^,(t) = ± V^, andz, = ^c + tGo(^c) = ±2 V^. 



That is, the singularity occurs precisely at the edge of the spec- 
trum. Furthermore, the resulting singularity is of the square 
root type. To see that, let us expand the characteristic equa- 
tion around the singular point. We get 



(13) 



It follows that, in the vicinity of the positive edge of the spec- 
trum z - Zc = '2. yfr, ^ - - ir'^"* Vz - Zc- Thus, as z moves 
towards Zc and is bigger than Zc, £, moves to on the real 
axis. When z becomes smaller than Zc, ^ moves away from 
along the imaginary axis. The imaginary part therefore ex- 
ists for z < Zc and yields a spectral density p{z) ~ yfzi^^^, 
in agreement with (fT2l i. This square root behaviour of the 
spectral density implies that in the vicinity of the edge of the 
spectrum, the number of eigenvalues in an interval of width 
s scales as Ns^^^, implying that the interlevel spacing goes as 

A^-2/3. 

In order to capture the fine structure of the level density 
in the vicinity of the edge, we need to take into account the 
1/A^ corrections. We wish to do that within the Dyson fluid 
picture. Trying to go beyond the large A^ limit at the level of 
Eq. (|9]l requires in particular handling the connected two-point 
function. We shall not do so directly, but shall rely on time 
dependent orthogonal polynomials, i.e., extend to our time- 
dependent setting tools that are familiar in random matrix the- 
ory. This will allow us to obtain simple equations that gen- 
eralize Eq. (|9]l (albeit not for the one-particle density), take 
the form of viscid Burgers equations, and are valid for any A^. 
This procedure allows, in principle, the calculation of all the 
correlation functions for finite A^ and f. 

For the considered ensemble, the relevant set of orthogonal 
polynomials is that of Hermite polynomials, defined as 



h,(x)^(-l)'e^^''^^e-''-^''\ 
These admit the following useful integral representation 



(14) 



ht(x) = i-iN)' 




dqq e 



)<. -NIlcf-VixqN 



(15) 



We shall also use the so-called monic polynomials, i.e. poly- 
nomials where the coefficient of highest order term is equal to 
unity: Tij^ix) = hk(x)/N'' - Y]'i=i(x " -^/)' with Xj denoting the 
(real) zeros of the Hermite polynomials. In the random walk 
described above, the probability distribution retains its form 
at all instants of time. It is then easy to construct polynomi- 
als that remain orthogonal with respect to the time-dependent 
measure exp(^-Nx^ /(2t)\: all that needs to be done is to re- 
place N hy N/t in Eq. ( fTsT l. One obtains then 



7Tk(x, t) 

which satisfy 




dqq e 



•J — Cx 



dxe. 2r jr^(x, T)nm(x, t) = 6„mcl, 



(16) 



(17) 
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with cl ^ n\ V2^T«+i/2^ The monic character of the 7r„'s is 
not affected by the time dependence. By using the integral 
representation (fT6] i. it is easy to show that the 7r„(jt:, t)'s satisfy 
the following equation 



dr7T„(x, t) = -V,d'^^7T„(x, t), 



(18) 



with v, - 1 /2N. This is like a diffusion equation with how- 
ever a negative diffusion constant (note however that n:„ is an 
analytic function of x, and that n„{-iy, r), with y real, satisfies 
a diffusion equation with a positive constant). At this point, 
one may perform an inverse Cole-Hopf transform, i.e. , define 
the new function 



fk(z, t) = 2vsd^ In TTkiz, r) 



1 * 

-y- 



1 



XiiT) 



(19) 



with Imz 0. The resulting equation for is the viscid Burg- 
ers Jsl] equation 



drfkiz, t) + Mz, T)djk(z, t) = -v.dtfkiz, t). 



(20) 



with V5 playing the role of a viscosity. 

The equation ( l20b is satisfied by all the functions f^. We 
shall focus now on the function associated to n^iz, r). One 
reason is that n^iz, t) is known to be equal to the average char- 
acteristic polynomial |8], i.e 



<det(z-//(r))) = ;rw(z,T), 



(21) 



and, in the large limit, ( 1 /A^)5; ln<det(z - H(t))} ^ 
(1/N)d,{\ndet{z-H{T))} = G(z). Thus /iv(z, t) coincides with 
the average resolvent G(z, r) in the large limit. In fact the 
structure of /ai, as clear from Eq. ( fT9l ), is very close to that of 
the resolvent, with its poles given by the zeros of the charac- 
teristic polynomial. Eq. ( l20l i for /iv(z, t) is exact. The initial 
condition, /iv(z, r = 0) = l/z, does not depend on A^, so that all 
the finite corrections are taken into account by the viscous 
term. Note however that this exact equation does not allow us 
to directly study the finite A^ corrections to the resolvent or the 
spectral density. 

We turn now to the study of the viscid Burgers equation for 
fN{z, t), i.e. Eq, ( |20l l for k - N,m the vicinity of the (moving) 
pre-shock, that is for x near Zc{t)- The Cole-Hopf transforma- 
tion used above provides a solution in term of n^ix, r) which 
allows us to study the effects of a small viscosity by perform- 
ing a saddle point approximation on the integral ( fTSI ). The 
saddle point equation coincides with the characteristic equa- 
tion discussed earlier (with the identification ^ — » -iq). How- 
ever this analysis breaks down in the vicinity of the shock: 
here the difficulties comes from the merging of the two sad- 
dle points associated with the square root singularity. A better 
analysis is then called for, and we shall rely in what follows 
on tools borrowed from the theory of turbulence 1 10]. Let us 
recall that in the vicinity of the edge of the spectrum, and in 
the inviscid limit. 



^ . ^ 1 1 r- 



(22) 



We set 



X = ZciT) + vf\s, Mx, t) = Zc(t) + vl'^XNis, t), (23) 

with Zc = drZc = ± 1 / V^- The particular scaling of the coor- 
dinate is motivated by the fact that near the square root sin- 
gularity the spacing between the eigenvalues scales as N^^^^. 
A simple calculation then yields the following equation for 
Xis, t) in the vicinity of Zc{t) - 2 V^"- 



1/3, 



drZc + V, 



dr ds 



,1/3 



(24) 



which, ignoring the small term of order v/", we can write as 

1 . dx 



d_ 



"2t3/2 ^ 2^ ^ ds 



or, setting - 25.5 In (^(i), as 



d^ip - a(s - SQ)(p - 0, 



= 0, 



4t3/2' 



(25) 



(26) 



with 5o(t) an arbitrary function of t. The general solu- 
tion of this equation is given by Airy functions fzT] 0,(.s) = 
e,A/(e,(5 - SQ)a^l^), where e,- for / = 0, 1, 2 are cubic roots of 1 
and Yjq - 0. In the case of the characteristic polynomial, 
the solution corresponds to eo = 1 ■ Finally, 



Aiia^l\s - so)) ■ 



(27) 



We complete our hydrodynamic description of the Dyson 
fluid by considering the Cauchy transforms of the monic or- 
thogonal polynomials 



Pkiz 



,r)=j- r 

2m J_o< 



dx 



TTkix, T)e 



(28) 



The motivation for doing so is that the average of the inverse 
characteristic polynomial is related to Pn-\ iflal : 



1 



det(z - H(t)) 



2m 



-pN-i(z, r). 



(29) 



where j is given after Eq. ( fTTb . To get the time dependent 
version of the pk's, we plug into Eq. ( |28] ) the integral repre- 
sentation of the time-dependent monic polynomials 7t^{x, t). 
A simple calculation then yields, for Imz > 0, 



Pkiz) = (-0' 



27rT Jo 



(30) 



and pki-z,T) - (-l)*^'pi(z, t). It can be shown by a direct 
calculation that the functions 



Pk(z,T) = e2r^'pi^(z,T) 



(31) 



satisfy the same recurrence relations and differential equations 
as the polynomials 7rt(z, t), in particular Eq. ( fTsT l. One may 
also verify that the function 



gkiz, t) = 2v,d, In Pkiz, t), 



(32) 
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analogous to fk(z,T) in Eq. ( fT9] ), satisfies a viscid Burgers 
equation: 

drgk + gkdzgk = -Vsdlgk- (33) 

Also this equation exhibits the phenomenon of universal pre- 
shock, but this time the scaling solution involves the two other 
solutions of the Airy equation that were mentioned earlier, 
namely (pi or 02, depending on the sign of the imaginary part 
of z, in agreement with LI 511 . 

We note that since nN(z) and Pn(z), and their first deriva- 
tives, are the building blocks of all relevant multipoint cor- 
relators (products of characteristic polynomials, products of 
inverse characteristic polynomials and mixed products iHgII ). 
all universal kernels originate implicitly from the dynamics of 
shocks in the viscid Dyson fluid. Also, the so-called "hard- 
edge" scaling (Bessel kernels) can be be viewed in the same 
way, as a pre-shock approaching a hard wall. We note also 
that other forms of universal shocks may appear, e.g., when 
the spectral density develops a cubic root singularity; then the 
universal character of the phenomenon is coded in a Pearcey 
function. This happens e.g. when two intervals supporting the 
hermitian spectrum merge as a function of an external param- 
eter 1 1 1, 12], or when two Airy type shock waves merge on a 
compact support 1^ Ull] . 

We expect the fluid dynamical picture presented here to 
hold also for other matrix-valued diffusion processes, in par- 
ticular for the multiplicative diffusions, corresponding to ran- 
dom walks involving products of complex matrices IItI [Tsll . 
Perhaps the best studied example is the random walk of 
unitary matrices [19], because of its relation with two di- 
mensional Yang-Mills theory with a large number of col- 



ors Nc f2d\, and the associated universality conjectured by 
Narayanan and Neuberger [13J . We have earlier emphasized 
the relevance of shock waves in this context i^, as provid- 
ing a simple mechanism for the Durhuus-Olesen f2T] order- 
disorder transition. Recently, Neuberger used the explicit ex- 
pressions of averages of the characteristic polynomial and its 
inverse, that were obtained through a character expansion, to 
prove that closely related functions satisfy Burgers equations 
with a spectral viscosity 1 /2Nc ll22i[23ll . The present analysis 
supports the general character of these results and may moti- 
vate further generalizations. In particular, we expect that pre- 
shock phenomena may be identified for non-compact random 
matrix models, e.g. in mesoscopic systems where the role of 
the time is played by the length of the wire ll24ll . Last but not 
least, this picture of shock formation finds analogies in other 
branches of physics, e.g., in the descritption of merging singu- 
larities in optics, leading to the so-called diffraction catastro- 
phes (see e.g. f2S]), or in the study of growth processes of the 
Kardar-Parisi-Zhang universality class and statistical proper- 
ties of the equilibrium shapes of crystals ll26ll . We believe 
that pushing Dyson's original concept of temporal dynamics 
for random matrices may allow for a better understanding of 
above analogies. 
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